#!/bin/tcsh -ef
####
#
#############################################################################
#                                                                           #
# Title: Generate Map from Saved Dataset                                    #
#                                                                           #
# 2dx.org, GNU Plublic License.                                             #
#                                                                           #
# Created..........: 02/20/2006                                             #
# Last Modification: 09/20/2006                                             #
# Author...........: 2dx.org                                                #
#                                                                           #
#############################################################################
#
# SORTORDER: 33
#
# MANUAL: This script makes use of MRC and CCP4 commands to generate the final projection map.
#
# DISPLAY: RESMIN
# DISPLAY: RESMAX
# DISPLAY: SYM
# DISPLAY: ALAT
# DISPLAY: tempfac
# DISPLAY: avrgamphsNUMBER
# DISPLAY: avrgamphsRESOL
# DISPLAY: zstarwin
# DISPLAY: zstarrange
# DISPLAY: Calc_from_zstarrange
# DISPLAY: zstarrange_real
# DISPLAY: realang
# DISPLAY: realcell
# DISPLAY: merge_modus
# DISPLAY: det_tilt
# DISPLAY: merge_modus
# DISPLAY: diffmap_doit
# DISPLAY: diffmap_source
# DISPLAY: diffmap_weight
# DISPLAY: merge_ref_num
# DISPLAY: merge_comment_1
# DISPLAY: merge_register_date_1
# DISPLAY: merge_comment_2
# DISPLAY: merge_register_date_2
# DISPLAY: merge_comment_3
# DISPLAY: merge_register_date_3
# DISPLAY: merge_comment_4
# DISPLAY: merge_register_date_4
# DISPLAY: merge_comment_5
# DISPLAY: merge_register_date_5
# DISPLAY: merge_comment_6
# DISPLAY: merge_register_date_6
# DISPLAY: merge_comment_7
# DISPLAY: merge_register_date_7
# DISPLAY: merge_comment_8
# DISPLAY: merge_register_date_8
# DISPLAY: merge_comment_9
# DISPLAY: merge_register_date_9
#
#$end_local_vars
#
set bin_2dx = ""
set proc_2dx = ""
#
set tempkeep = ""
set RESMIN = ""
set RESMAX = ""
set CS = ""
set KV = ""
set ALAT = ""
set realang = ""
set realcell = ""
set rot180 = ""
set revhk = ""
set rot90 = ""
set beamtilt = ""
set zstarwin = ""
set zstarrange = ""
set Calc_from_zstarrange = ""
set zstarrange_real = ""
set tempfac = ""
set SYM = ""
set avrgamphsNUMBER = ""
set avrgamphsRESOL = ""
set merge_modus = ""
set det_tilt = ""
set merge_modus = ""
set diffmap_doit = ""
set diffmap_source = ""
set diffmap_weight = ""
set merge_ref_num = ""
set merge_comment_1 = ""
set merge_comment_2 = ""
set merge_comment_3 = ""
set merge_comment_4 = ""
set merge_comment_5 = ""
set merge_comment_6 = ""
set merge_comment_7 = ""
set merge_comment_8 = ""
set merge_comment_9 = ""
set merge_register_date_1 = ""
set merge_register_date_2 = ""
set merge_register_date_3 = ""
set merge_register_date_4 = ""
set merge_register_date_5 = ""
set merge_register_date_6 = ""
set merge_register_date_7 = ""
set merge_register_date_8 = ""
set merge_register_date_9 = ""
#
#$end_vars
#
echo "<<@progress: 1>>"
#
set scriptname = 2dx_genMergeMapRegister
\rm -f LOGS/${scriptname}.results
#
set ccp4_setup = 'y'
source ${proc_2dx}/initialize
#
set IS_2DX = yes
source ${proc_2dx}/2dx_makedirs 
#
set imagenumber = 1001
set imagename = "merge"
#
echo imagename = ${imagename}
echo imagenumber = ${imagenumber}
echo tempkeep = ${tempkeep}
echo RESMIN = ${RESMIN}
echo RESMAX = ${RESMAX}
echo CS = ${CS}
echo KV = ${KV}
echo ALAT = ${ALAT}
echo realang = ${realang}
echo realcell = ${realcell}
echo rot180 = ${rot180}
echo rot90 = ${rot90}
echo beamtilt = ${beamtilt}
echo zstarwin = ${zstarwin}
echo zstarrange = ${zstarrange}
echo tempfac = ${tempfac}
echo SYM = ${SYM}
echo avrgamphsNUMBER = ${avrgamphsNUMBER}
echo avrgamphsRESOL = ${avrgamphsRESOL}
echo diffmap_doit = ${diffmap_doit}
echo diffmap_source = ${diffmap_source}
echo diffmap_weight = ${diffmap_weight}
#
set ABANG = `echo $realang | awk '{s=180-$1} END {print s}'`
echo ABANG = ${ABANG}
set IAQP2 = 0
set IVERBOSE = 1
set phastepnum = 1
set phastep = 0.1
#
if ( ${ALAT} == "0" || ${ALAT} == "0.0" ) then
  ${proc_2dx}/protest "ALAT is not defined."
endif
set ALATnew = `echo ${ALAT} | awk '{ if ( $1 < 0 ) { s = -$1 } else { s = $1 }} END { print s }'`
if ( ${ALAT} != ${ALATnew} ) then
  set ALAT = ${ALATnew}
  echo "set ALAT = ${ALAT}" >> LOGS/${scriptname}.results
endif
#
set zwin = `echo ${ALAT} | awk '{ s = 1.0 / ( 2.0 * $1 ) } END { print s }'`
echo zwin = $zwin
#
if ( ${zstarwin} == "0" ) then
  set zstarwin = ${zwin}
  echo "set zstarwin = ${zstarwin}" >> LOGS/${scriptname}.results
endif
#
if ( ${Calc_from_zstarrange} == "y" ) then
  set zstarrange_real = `echo ${zstarrange} ${ALAT} | awk '{ s = 1.0 / ( $1 ) } END { print s }'`
  echo "set zstarrange_real = ${zstarrange_real}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "Calculating vertical resolution as ${zstarrange_real} Angstroems."
else
  set zstarrange = `echo ${zstarrange_real} ${ALAT} | awk '{ s = 1.0 / ( $1 ) } END { print s }'`
  echo "set zstarrange = ${zstarrange}" >> LOGS/${scriptname}.results
  ${proc_2dx}/linblock "Calculating zstarrange as ${zstarrange} (with 0.5 = Nyquist resolution)."
endif
#
set zmin = `echo ${zstarrange} | awk '{s = -$1} END {print s}'`
set zminmax = `echo ${zmin},${zstarrange}`
#
if ( ${merge_modus} == "2D" ) then
  set zminmax = "-0.01,0.01"
else
  set zminmax = "-0.5,0.5"
endif
#
echo zminmax = ${zminmax}
#
#???? Check this ???
set SCL = 1
echo "SCL = ${SCL}"
#
# contrast for grey plot
set scale = 1
echo "scale = ${scale}"
#
if ( ${merge_ref_num} != "0" ) then
  echo ":: "
  ${proc_2dx}/linblock "Using merged dataset from register ${merge_ref_num}."
  echo ":: "
  if ( ${merge_ref_num} == "1" ) then
    ${proc_2dx}/linblock "Comment: ${merge_comment_1}."
    ${proc_2dx}/linblock "Date: ${merge_register_date_1}."
  endif
  if ( ${merge_ref_num} == "2" ) then
    ${proc_2dx}/linblock "Comment: ${merge_comment_2}."
    ${proc_2dx}/linblock "Date: ${merge_register_date_2}."
  endif
  if ( ${merge_ref_num} == "3" ) then
    ${proc_2dx}/linblock "Comment: ${merge_comment_3}."
    ${proc_2dx}/linblock "Date: ${merge_register_date_3}."
  endif
  if ( ${merge_ref_num} == "4" ) then
    ${proc_2dx}/linblock "Comment: ${merge_comment_4}."
    ${proc_2dx}/linblock "Date: ${merge_register_date_4}."
  endif
  if ( ${merge_ref_num} == "5" ) then
    ${proc_2dx}/linblock "Comment: ${merge_comment_5}."
    ${proc_2dx}/linblock "Date: ${merge_register_date_5}."
  endif
  if ( ${merge_ref_num} == "6" ) then
    ${proc_2dx}/linblock "Comment: ${merge_comment_6}."
    ${proc_2dx}/linblock "Date: ${merge_register_date_6}."
  endif
  if ( ${merge_ref_num} == "7" ) then
    ${proc_2dx}/linblock "Comment: ${merge_comment_7}."
    ${proc_2dx}/linblock "Date: ${merge_register_date_7}."
  endif
  if ( ${merge_ref_num} == "8" ) then
    ${proc_2dx}/linblock "Comment: ${merge_comment_8}."
    ${proc_2dx}/linblock "Date: ${merge_register_date_8}."
  endif
  if ( ${merge_ref_num} == "9" ) then
    ${proc_2dx}/linblock "Comment: ${merge_comment_9}."
    ${proc_2dx}/linblock "Date: ${merge_register_date_9}."
  endif
  if ( -d REGISTERS/Reg_${merge_ref_num} ) then
    if ( ${merge_modus} == "2D" ) then
      if ( ! -e REGISTERS/Reg_${merge_ref_num}/merge2D_MRClefthanded.mtz ) then
        ${proc_2dx}/protest "ERROR: merge2D_MRClefthanded.mtz not existing in register ${merge_ref_num}."
      endif
      set infile = REGISTERS/Reg_${merge_ref_num}/merge2D_MRClefthanded.mtz
    else
      if ( ! -e REGISTERS/Reg_${merge_ref_num}/merge3D_MRClefthanded.mtz ) then
        ${proc_2dx}/protest "ERROR: merge3D_MRClefthanded.mtz not existing in register ${merge_ref_num}."
      endif
      set infile = REGISTERS/Reg_${merge_ref_num}/merge3D_MRClefthanded.mtz
    endif
  else
    ${proc_2dx}/protest "ERROR: Register ${merge_ref_num} does not contain data."
  endif
  echo ":: "
  #
else
  ${proc_2dx}/linblock "Using the most recent merged dataset (from register 0)."
  if ( ${merge_modus} == "2D" ) then
    set infile = merge2D_MRClefthanded.mtz
  else
    set infile = merge3D_MRClefthanded.mtz
  endif
endif
#
echo "# IMAGE: ${infile} <MTZ: input data>"  >> LOGS/${scriptname}.results
#
#############################################################################
${proc_2dx}/linblock "sourcing sym2spsgrp_sub.com"
#############################################################################
#
source ${proc_2dx}/2dx_sym2spcgrp_sub.com
#
echo SYM = ${SYM}
echo spcgrp = ${spcgrp}
echo CCP4_SYM = ${CCP4_SYM}
#
#############################################################################
${proc_2dx}/linblock "mtzdump - to look at ${imagename}_MRClefthanded.mtz"
#############################################################################
#
${bin_ccp4}/mtzdump hklin ${infile} << eot
NREF 100
END
eot
#
#############################################################################
${proc_2dx}/linblock "fft - to transform ${infile} into SCRATCH/scratch1.map"
#############################################################################
#
echo "<<@progress: 20>>"
#
# contrast for grey plot
set scale = 1
echo "scale = ${scale}"
#
\rm -f SCRATCH/scratch1.map
#
if ( ${merge_modus} == "3D" ) then
  #############################################################################
  #############################################################################
  #############################################################################
  # Here for 3D volume generation:
  #############################################################################
  #############################################################################
  #############################################################################
  #
  if ( ${CCP4_SYM} == "P2" ) then
    #
    if ( ${revhk} == 'n' ) then
      set AXIS = "Z,X,Y"
    else
      ${proc_2dx}/linblock "four: revhk = ${revhk}"
      set AXIS = "X,Z,Y"
    endif 
    #
    ${bin_ccp4}/fft hklin ${infile} mapout SCRATCH/scratch1.map  << eot
LABIN F1=F  PHI=PHI ##
AXIS ${AXIS}
SCALE F1 ${scale} ${tempfac}
SYMMETRY P2
RESOLUTION ${RESMIN} ${RESMAX}
TITLE ${imagename}, ${SYM}, res=${RESMAX}, T=${tempfac}, ${date}
GRID 200 200 200
XYZLIM 0 199 0 199 0 199
RHOLIM 250.0
HKLMAX 50 50 50
END
eot
    #
  else
    #
    if ( ${revhk} == 'n' ) then
      set AXIS = "X,Y,Z"
    else
      ${proc_2dx}/linblock "four: revhk = ${revhk}"
      set AXIS = "Y,X,Z"
    endif
    #
    ${bin_ccp4}/fft hklin ${infile} mapout SCRATCH/scratch1.map  << eot
LABIN F1=F  PHI=PHI ##
AXIS ${AXIS}
SCALE F1 ${scale} ${tempfac}
SYMMETRY ${CCP4_SYM}
RESOLUTION ${RESMIN} ${RESMAX}
TITLE ${imagename}, ${SYM}, res=${RESMAX}, T=${tempfac}, ${date}
GRID 200 200 200
XYZLIM 0 199 0 199 0 199
RHOLIM 250.0
HKLMAX 50 50 50
END
eot
    #
  endif
  #
  echo "<<@progress: 50>>"
  #
  echo "# IMAGE: SCRATCH/scratch1.map <MAP: temporary volume 1>" >> LOGS/${scriptname}.results
  #
  #############################################################################
  ${proc_2dx}/linblock "maprot - to rescale map into correct dimensions"
  #############################################################################
  #
  set cellx = `echo ${realcell} | cut -d\, -f1`
  set cellxm1 = `echo ${cellx} | awk '{s=$1-1} END {print s}'`
  echo ":cellx = ${cellx},    cellxm1 = ${cellxm1}"
  #
  set celly = `echo ${realcell} | cut -d\, -f2`
  set cellym1 = `echo ${celly} | awk '{s=$1-1} END {print s}'`
  echo ":celly = ${celly},    cellym1 = ${cellym1}"
  #
  set ALATm1 = `echo ${ALAT} | awk '{s=$1-1} END {print s}'`
  echo ":ALAT = ${ALAT},      ALATm1 = ${ALATm1}"
  #
  \rm -f SCRATCH/scratch2.map
  ${bin_ccp4}/maprot mapin SCRATCH/scratch1.map wrkout SCRATCH/scratch2.map << eot
MODE FROM
CELL WORK ${realcell} ${ALAT} 90.0 90.0 ${realang}
GRID WORK ${cellx} ${celly} ${ALAT}
XYZLIM 0 ${cellxm1} 0 ${cellym1} 0 ${ALATm1}
SYMM WORK ${CCP4_SYM}
AVER
OMAT
  1.000 0.000 0.000
  0.000 1.000 0.000
  0.000 0.000 1.000
  0.000 0.000 0.000
eot
  #
  echo "<<@progress: 80>>"
  #
  echo "# IMAGE: SCRATCH/scratch2.map <MAP: temporary volume 2>" >> LOGS/${scriptname}.results
  #
  #############################################################################
  ${proc_2dx}/linblock "mapmask - to set units to 2x2x1"
  #############################################################################
  #
  \rm -f volume.map
  ${bin_ccp4}/mapmask mapin SCRATCH/scratch1.map mapout volume.map << eof
AXIS X,Y,Z
scale factor 1
xyzlim -1.0 1.0 -1.0 1.0 -0.5 0.5
END
eof
  #
  echo "# IMAGE-IMPORTANT: volume.map <MAP: Final Volume>" >> LOGS/${scriptname}.results
  #
  echo "<<@progress: 100>>"
  #
  exit
  #
else
  #
  #############################################################################
  #############################################################################
  #############################################################################
  # Here for 2D work:
  #############################################################################
  #############################################################################
  #############################################################################
  #
  if ( ${CCP4_SYM} == "P2" ) then
    #
    if ( ${revhk} == 'n' ) then
      set AXIS = "Z,X,Y"
    else
      ${proc_2dx}/linblock "four: revhk = ${revhk}"
      set AXIS = "X,Z,Y"
    endif 
    #
    ${bin_ccp4}/fft hklin ${infile} mapout SCRATCH/scratch1.map  << eot
LABIN F1=F  PHI=PHI ##
PROJECTION
AXIS ${AXIS}
SCALE F1 ${scale} ${tempfac}
SYMMETRY P2
RESOLUTION ${RESMIN} ${RESMAX}
TITLE ${imagename}, ${SYM}, res=${RESMAX}, T=${tempfac}, ${date}
GRID 200 2 200
XYZLIM 0 199 0 0 0 199
RHOLIM 250.0
HKLMAX 50 50 50
END
eot
    #
  else
    #
    if ( ${revhk} == 'n' ) then
      set AXIS = "X,Y,Z"
    else
      ${proc_2dx}/linblock "four: revhk = ${revhk}"
      set AXIS = "Y,X,Z"
    endif
    #
    ${bin_ccp4}/fft hklin ${infile} mapout SCRATCH/scratch1.map  << eot
LABIN F1=F  PHI=PHI ##
PROJECTION
AXIS ${AXIS}
SCALE F1 ${scale} ${tempfac}
SYMMETRY ${CCP4_SYM}
RESOLUTION ${RESMIN} ${RESMAX}
TITLE ${imagename}, ${SYM}, res=${RESMAX}, T=${tempfac}, ${date}
GRID 200 200 20
XYZLIM 0 199 0 199 0 0
RHOLIM 250.0
HKLMAX 50 50 50
END
eot
    #
  endif
endif
#
#############################################################################
# The remainder is also for 2D projection maps only:
#############################################################################
#
echo "<<@progress: 50>>"
#
#############################################################################
${proc_2dx}/linblock "extend - to extend SCRATCH/scratch1.map into ${imagename}-${SYM}.rec.mrc"
#############################################################################
#
\rm -f SCRATCH/${imagename}-${SYM}.map
#
if ( ${CCP4_SYM} == "P2" ) then
  #
  \rm -f SCRATCH/TMP001.map
  #
${bin_ccp4}/extends mapin SCRATCH/scratch1.map mapout SCRATCH/TMP001.map << eof
XYZLIM 0 399 0 0 0 399
KEEP
END
eof
  #
  \rm -f SCRATCH/TMP001.mrc
  #
  ${bin_2dx}/labelh.exe << eot
SCRATCH/TMP001.map
13
SCRATCH/${imagename}-${SYM}.map
eot
  #
else
  #
${bin_ccp4}/extends mapin SCRATCH/scratch1.map mapout SCRATCH/${imagename}-${SYM}.map << eof
XYZLIM 0 399 0 399 0 0
KEEP
END
eof
#
endif
#
echo "# IMAGE: SCRATCH/${imagename}-${SYM}.map <temporary map>" >> LOGS/${scriptname}.results
#
echo "<<@progress: 55>>"
#
#############################################################################
${proc_2dx}/linblock "npo - to create a line plot ${imagename}-${SYM}.plt"
#############################################################################
#
\rm -f ${imagename}-${SYM}.plt
#
${bin_ccp4}/npo  MAPIN  SCRATCH/${imagename}-${SYM}.map  PLOT  ${imagename}-${SYM}.plt  << eof
TITLE Merged Result
MAP SCALE 0.4
CONTRS -500 TO 500 BY 12
MODE BELOW -40 DASHED 1 0.15 0
SECTS 0 0
PLOT
END
eof
#
echo "<<@progress: 85>>"
#
#############################################################################
${proc_2dx}/linblock "laserplot - to create PS/${imagename}MAP-${SYM}.ps"
#############################################################################
#
\rm -f PS/${imagename}MAP-${SYM}.ps
${bin_2dx}/laserplot.exe -outputfile=PS/${imagename}MAP-${SYM}.ps ${imagename}-${SYM}.plt
#
\rm -f ${imagename}-${SYM}.plt
#
echo "# IMAGE: PS/${imagename}MAP-${SYM}.ps <PS: ${SYM}-symmetrized final map plot>" >> LOGS/${scriptname}.results
#
echo "<<@progress: 90>>"
#
#############################################################################
${proc_2dx}/linblock "LABEL - to create a clean MRC file format instead of CCP4"
#############################################################################
#
\rm -f ${imagename}-${SYM}.mrc
#
${bin_2dx}/labelh.exe << eot
SCRATCH/${imagename}-${SYM}.map 
2
${imagename}-${SYM}.mrc
1,0
0
eot
#
echo "# IMAGE: ${imagename}-${SYM}.mrc <${SYM}-symmetrized final map>" >> LOGS/${scriptname}.results
#
echo "<<@progress: 100>>"
#
#############################################################################
#                                                                           #
${proc_2dx}/linblock "Done."
#                                                                           #
#############################################################################
#
